Quantum phase transitions in the Kane-Mele-Hubbard model 
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We study the two-dimensional Kane-Mele-Hubbard model at half filling by means of quantum 
Monte Carlo simulations. We present a refined phase boundary for the quantum spin liquid. The 
topological insulator at finite Hubbard interaction strength is adiabatically connected to the ground- 
state of the Kane-Mele model. In the presence of spin-orbit coupling, magnetic order at large Hub- 
bard U is restricted to the transverse direction. The transition from the topological band insulator 
to the antiferromagnetic Mott insulator is in the universality class of the three-dimensional XY 
model. The numerical data suggest that the spin liquid to topological insulator and spin liquid to 
Mott insulator transitions are both continuous. 
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I. INTRODUCTION 

Topological insulators have attracted significant atten- 
tion in recent years? 1 , especially since their experimen- 
tal realization^ Whereas the existence of the topological 
state and many of its consequences can be understood in 
terms of exactly solvable, noninteracting models, the in- 
terplay of a topological band structure and electronic cor- 
relations has become a very active field of research. The 
corresponding interacting models do not have general ex- 
act solutions, which has made computational methods 
one of the most important tools. A possible experimen- 
tal route to the strongly correlated regime is based on 
optical lattices^ 

The Z2 topological band insulator (TBI), or quan- 
tum spin-Hall insulator, closely related to the integer 
quantum Hall effect^ can be realized in the Kane-Mele 
(KM) modeli^£ The latter describes electrons (or Dirac 
fermions) on the two-dimensional (2D) honeycomb lat- 
tice, with nearest-neighbor hopping and spin-orbit cou- 
pling. Originally motivated by graphene^ the spin-orbit 
coupling turned out to be much too small in this mate- 
rial for topological effects to be observable. However, the 
KM model and its extension, the Kane-Mele-Hubbard 
(KMH) model turn out to be a very useful theoretical 
framework. In particular, the honeycomb lattice geome- 
try provides a direct connection to the recently discovered 
quantum spin liquid (QSL) phase of the Hubbard model 
on the same latticed For the latter, the Dirac spectrum 
with vanishing density of states at the Fermi level leads 
to a Mott transition at a finite critical Hubbard U, and 
the QSL phase lies between a semimetal and a magnetic 
insulator** Finally, the symmetries of the KMH model 
permit the application of powerful quantum Monte Carlo 
(QMC) methods without a sign problem,— £ so that exact 
results can be obtained. 

The phase diagram of the KMH model has been de- 
rived from QMC simulations^^ and numerical results for 
the extent of the QSL were presented in Ref. 0- At any 



nonzero spin-orbit coupling, the semimetal is replaced by 
the Z2 TBI. In contrast, the gapped QSL is found to be 
stable up to a finite critical value of the spin-orbit inter- 
action. Finally, the magnetic transition of the Hubbard 
model, between the QSL and an antiferromagnetic Mott 
insulator (AFMI), is supplemented with a similar transi- 
tion between the TBI and the AFMI in a potentially dif- 
ferent universality class. On a qualitative level, certain 
aspects of the phase diagram were obtained for example 
in mean-field theory^ as well as with cluster method a 10 ' 11 
and variational QMC— 

The understanding of the KMH model is not com- 
plete. Many of the open questions are related to the 
perhaps most intriguing aspect of the model, namely the 
QSL phase. The recent results from approximate cluster 
methods for parameters in the QSL region of the exact 
phase diagram inaccurately suggest a rather complete un- 
derstanding of this exotic phase. However, strictly speak- 
ing, any cluster method breaks translational symmetry, 
so that a true QSL phase is excluded from the outset. In 
this light, conclusions such as the absence of edge states, 
or the closing of the single-particle gap across the tran- 
sition to the TBI are not surprising, as the QSL phase is 
replaced in these studies by a simple band insulator (a va- 
lence bond crystal) . The large correlation lengths (small 
gaps) observed in the QSL phase in the Hubbard model^ 
highlight the necessity of careful interpretation of the re- 
sults obtained by cluster approximations in the context 
of the QSL. Interesting connections between TBIs and 
QSLs are discussed in Ref. [H. 

The purpose of this paper is threefold. First, we 
present a much more detailed account of the QMC cal- 
culations underlying the phase diagram shown in Ref. @- 
Second, we extend the number of points in parameter 
space and the observables calculated, in order to pro- 
vide additional insight. We also present a refined phase 
boundary for the QSL phase. Third, we use the QMC 
method to investigate the quantum phase transitions, es- 
pecially in the light of recent theoretical prcdictions i 14 i 15 
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We show that the TBI-AFMI transition is in the ex- 
pected 3D XY universality class, and provide evidence 
for the continuous nature of the QSL-TBI and the QSL 
AFMI quantum phase transitions. In contrast to earlier 
workji we only consider bulk properties. We also pro- 
vide an overview of recent work on correlation effects in 
topological insulators with a focus on the KMH model. 

The paper is organized as follows. In Scc.[H]wc briefly 
review the model. Details about the QMC method are 
presented in Sec. IIIII Section IIVI contains our numeri- 
cal results, beginning with the refined phase diagram, 
and followed by a detailed account of the various quan- 
tum phase transitions. We end with conclusions and an 
overview of open questions in Sec. |V] 



II. MODEL 

The Hamiltonian of the KMH model can be written in 
the form H = Hkm + Hjj, where 

Hkm = -* c i c J +iX J2 v iA aZ °o ' 



U 



^ = yE( c k-i) 2 



(1) 



Here c] = (cj ^, c\ jj is a spinor of electron creation oper- 
ators, i is the position of a lattice site on the honeycomb 
lattice, denotes a pair of nearest neighbors, and 

is a pair of next-nearest-neighbor lattice sites; a z 
is a Pauli matrix, and u%j = ±1 depending on whether the 
hopping path defined by the nearest-neighbor bonds con- 
necting sites i and j bends to the right or to the left. The 
complex next-nearest-neighbor hopping term in Eq. ([T]) 
can be related to the spin-orbit interaction in graphene 
and accounts for a spin-dependent staggered magnetic 
field.— The choice of writing the interaction term Hjj in 
an SU(2) invariant form is related to previous work on 
the Hubbard model on the honeycomb lattice^ in which 
it was essential to build this symmetry into the QMC 
method, see also Sec. IIIII In the presence of the spin- 
orbit term, the SU(2) spin rotation symmetry is reduced 
to a U(l) symmetry Throughout this work, we use peri- 
odic boundary conditions so that there are no edges, and 
take t as the unit of energy. The number of unit cells 
in each direction is denoted by L, the total number of 
unit cells is L 2 , and the total number of lattice sites is 
N = 2L 2 . The lattice sizes used satisfy L = 31 with I 
integer, and range from L = 3toL=18. We exclusively 
consider the case of a half-filled band. 

A possible Rashba term is neglected from the outset, 
because it would cause a sign problem in the QMC simu- 
lations. However, a small but finite Rashba coupling does 
not destroy the TBI state of the KM model^ The nonin- 
teracting case U — has been solved in the original paper 
by Kane and Mele^ Most importantly, the groundstate 
is a Z2 TBI for any finite spin-orbit coupling A. The KM 
model is closely related to a spinlcss model proposed by 



Haldanc which shows a quantum Hall effect and breaks 
time reversal invariance (TRI)J^ Combining two copies 
of the Haldane model gives the KM model exhibiting the 
quantum spin-Hall effect and preserving TRIi^i^ 

Hamiltonian ([T]) has been studied by means of 
mean-field and analytical approachesj2iiZr-2£ QMC 
simulations jZiSiiS the variational cluster approach^ clus- 
ter dynamical mean-field theory^ and field theory 14 i 15 
A more detailed discussion of previous results will be 
given in Sec. IIVI 



III. METHOD 

We employ a projective auxiliary-field determinant 
QMC algorithm similar to Ref. [H, which has previously 
been applied to the KMH modeli^ The method is based 
on the relation 

(..PI..)- to < V| e % T) ' T> 0) 

for the expectation value of an operator O, with a trial 
wave function | \&t) that is required to be nonorthogo- 
nal to the groundstate |^o)- It is beyond the scope of 
this article to describe the details of the algorithm, and 
the interested reader is referred to Ref. [2f|. Instead, we 
concentrate on aspects specific to the calculations pre- 
sented here, namely the choice of the trial wave function, 
the Hubbard-Stratonovich (HS) transformation, and the 
absence of the minus-sign problem for a half-filled band. 



A. Trial wave function 

In order to simplify the implementation, the trial wave 
function is taken to be a single Slater determinant, and 
can hence always be written in terms of the groundstate 
of a single-particle Hamiltonian Hi. There are many 
possible choices for |^t)- One can for example decide to 
optimize the overlap with the groundstate at the expense 
of symmetries ^ Here we have preserved symmetries, and 
have chosen |\&t) to be the groundstate of the KM model, 
which is defined by the first line of Eq. ([TJ. For A 7^ 0, the 
groundstate of the nonintcracting problem at half filling 
is insulating. Hence, the trial wave function is nondegen- 
erate and has all the symmetries of the Hamiltonian. At 
A = 0, the situation is more delicate. For the considered 
lattice sizes, L = 3/, the two nonequivalent Dirac points 
are located at the Fermi surface, and the groundstate of 
the noninteracting model at half filling is four-fold de- 
generate in each spin sector. We lift this degeneracy by 
means of a twist in the boundary condition in Ht in the 
direction a\ = (1,0), 



Ht = ~t ctcj exp 

(i,3) 



2ni 
5>n 



A-dl 



(3) 
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with A = Qai/L. The twist preserves translation sym- 
metry, so that the total momentum remains a good quan- 
tum number. In particular, for an infinitesimal twist, the 
groundstate has vanishing total momentum. Because fi- 
nite values of $ lead to a breaking of the C3 lattice sym- 
metry, the trial wave function cannot be classified ac- 
cording to the irreducible representation of this group at 
A = 0. After lifting possible degeneracies, |^t) corre- 
sponds to the nondegencratc groundstate of H^. This 
implies the relation 

lim e - 0(H ^- E ^ = |* t )(*t| , (4) 

O— >-oo 

where £r is the corresponding groundstate energy. 

B. Hubbard-Stratonovich transformation 

We choose a HS transformation of the Hubbard term 
Hjj that couples to the total density m = + 71,4., 
thereby conserving the SU(2) spin symmetry for every 
field configuration. In principle, HS transformations that 
couple to the z-component of spin are also possible. How- 
ever, at low temperatures, it is often difficult to restore 
the A = SU(2) spin symmetry of the total Hamilto- 
nian by stochastic sampling. An SU(2) spin s ym metric 
transformation was previously used in Refs. I23U24 

After a Trotter decomposition with imaginary time 
step At to isolate the interaction term, our HS trans- 
formation for a general operator O reads 

e -Ar0 2 = J- ^(ly^v(l)O +0 (At 4 ), (5) 
l=±l.±2 

I 



with the two functions and rj(l) of the auxiliary field 
I (with I = ±1, ±2) taking on the values 

7 (±l) = (l + V6/3)/4, 77(±1) = ±^2(3- V6), 

7(±2) = (1 - \/6/3)/4 , r?(±2) = ±^2(3 + ^6). 

(6) 

Equation ([5]) is an approximation to the Gaussian inte- 
gral and introduces an overall systematic error of the or- 
der At 3 , which is negligible in comparison to the Trotter 
error of order At 2 . The major advantage of this approx- 
imation is that we can avoid using continuous auxiliary 
fields while retaining spin rotation symmetry. For the 
Hubbard interaction, we have O = V / Uj2(n t +714.-1) 
and it is understood that the HS fields acquire space and 
time indices, 1 1— > li, T . 



C. Absence of a sign problem 

We prove the absence of the minus-sign problem for the 
projective QMC method at half filling. With the Trot- 
ter decomposition, choice of trial wave function and HS 
transformation, the denominator of Eq. ((2J factors into 
spin-up and spin-down determinants, 



<* T in e ~ Ar 



-AtH km -AtH l 



|*t) = Tr 



lim e ~ e ^- E ^ TT e -^KM e -Arff L 
e-s-oo 



!«« Ennn^R 



{l iyT } a t=1 i 



. (7) 

where we have used Eq. (j4|) to introduce a trace, {h,r} denotes an auxiliary-field configuration, and with the weights 



W a = Tr [e 9 ^ cxp [ - 9 £ c^W^cA fj exp [ - At £ c\ a (A t + cA 

L L ij > r=l L ij ' 

x cxp r)(li, T )(n i<7 - 1/2) 



(8) 



Here we introduced the notation H T = J2ija c lA h T(®)]ij c ja and H KM = J2ija c iA A t + A^.^^c^. Proving the 
absence of a negative sign problem at half filling amounts to showing that WZ = W\.. Since the trace in Eq. ([5J is 
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over one spin sector, we drop the spin index on the fermion operators to lighten the notation and obtain 



W* = Tr 



e° E - exp [ - 6 £ c\ [h* T ($)] i3 .c : . ] f[ exp { - At £ c t (A t * + A^c . | 

x expj - iy/ArU/2 V(kr)(ni - 1/2) 



Tr 



T=l 



e e£T exp j - e J2 ^T^wMr+M- r II ex p { - At E c i( A * + ^.tJijC-i)* 4 -^} 

exp i iy/ArU/2 ^ J?(Ii,r)(l - «i - 1/2 



(9) 



The second line follows from the canonical transforma- 
tion C; — > (— l) l ct, where the phase factor (—1)* takes 
the value 1 (—1) on sublattice A (B). The Hamiltonian 
i?T which generates the trial wave function has nonva- 
nishing matrix elements only between sites on opposite 
sublattices. Hence (— = —1 and 



Ci[>*(*)U-i) <+i <4 = cta*)]^ = ci.[/i T ($)],. iC . 



Similarly, for the hopping term, 



(10) 



(11) 



Since the spin-orbit term involves hopping between sites 
on the same sublattice, we have 



c^a.tM-1) 



i+3 c \ 



-4 (Ax,t) 



ji C i 



4 ( A Kl)u c i ■ 



Using Eqs. (fT0 l) -(fT2" j) . one sees that indeed 



Wf = W_ 



(12) 



(13) 



so that no sign problem exists at the particle-hole sym- 
metric point of the KMH model in the present formula- 
tion of the QMC algorithm^ The underlying reason is 
time reversal symmetry, which implies A\ t i = —A\ t f. 



D. Measurements 

For a given auxiliary-field configuration, we have to 
solve a free-electron Hamiltonian with external fields that 
vary in time and space. Consequently, Wick's theorem 
holds, and it is sufficient to compute the single-particle 
Green functions 

G a (i,3,r,T') = - (*o|Tc iiCT (r)4 iCT (r') |*„) (14) 

to calculate arbitrary correlation functions. For the cal- 
culation of G a we have followed Ref. The single- 
particle gap A sp at the Dirac point and the (staggered) 
spin gap A s at q = are extracted from fits to the cor- 
responding Green functions^ 



E. Projection parameter and Trotter discretization 

The projective algorithm involves two numerical pa- 
rameters, namely the projection parameter 9 and the 
Trotter time step At. Both parameters were chosen such 
that their influence on the results is smaller than the sta- 
tistical errors. Explicitly, we used Art — 0.05 or 0.1, and 
9t = 40-60. 



IV. RESULTS 

In order to better orient the discussion, we first present 
the phase diagram of the KMH model in Sec. IIV Al to- 
gether with a review of recent work, before elaborating 
on the various quantum phase transitions. 



A. Phase diagram 

Figure [T] shows the groundstate phase diagram of the 
KMH model, as obtained from QMC simulations. In ad- 
dition to the three phases of the Hubbard model on the 
honeycomb lattice, the spin-orbit coupling introduces a 
Z2 TBI. The gapless SM phase exists only at A = 0. 
Whereas the TBI and the QSL phase are fully gapped 
(finite single-particle gap A sp and spin gap A s ), the mag- 
netic phase has A sp > but A s = 0. Here all gaps refer 
to the bulk, and are not to be confused with the metallic, 
gapless edge states of the TBI phase of the KMH model. 

To the best of our knowledge, the QSL phase is char- 
acterized by the absence of any local order parameter 
which would reflect a broken-symmetry state. It can 
hence be regarded as a genuine Mott insulating state, 
which should be stable with respect to small perturba- 
tions such as spin-orbit coupling. In the case of an odd 
number of electrons per unit cell, the generalization of the 
Lieb-Schultz-Mattis theorem to two dimensions^ sug- 
gests the presence of topological order in the most gen- 
eral sense. Since the half-filled honeycomb lattice has two 
electrons per unit cell, this topological ordering still has 
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FIG. 1: (Color online) Groundstate phase diagram of the 
Kane-Mele-Hubbard model as obtained from QMC simula- 
tions. The four phases are a Z2 topological band insula- 
tor (TBI) with nonzero single-particle (spin) gap A sp > 
(A s > 0), a semimetal (SM, A sp = A s = 0) existing at A = 0, 
a quantum spin liquid (QSL, A sp > 0, A s > 0), and an an- 
tiferromagnetic Mott insulator (AFMI, A sp > 0, A s = 0). 
Magnetic order in the z direction exists in the AFMI at A = 0, 
but can be excluded for X/t > 0.002 and all values of U/t 
shown. Lines are quadratic fits to the QMC data points. 



to be numerically demonstrated or refuted. The underly- 
ing SU(2) x SU{2)/Z 2 symmetry of the Hubbard model 
on the honeycomb lattice has led to the prediction of a 
Z-2 x Z2 QSL with mutual spin-charge statistics.— Sub- 
lattice pairing states have been put forward by various 
authors to account for the QSL phasei 28 ' 29 A canonical 
consequence of the above topologically ordered phases 
is that, assuming a continuous phase transition, the 
magnetically ordered phase would not be a simple Ncel 
state, thus leading to conjectures that can be tested 
numericallyj 28 ' 29 Finally, in the presence of spin-orbit 
coupling, the possibility of the emergence of a topologi- 
cal Mott insulating phase, in which the spinons carry the 
topological character of the phase, remains i 30 i 31 For the 
KMH model, recent theoretical suggestions include a Z2 
QSLi£ and a chiral QSL.^S 

The boundary of the magnetic phase is obtained from 
the onset of long-range antifcrromagnctic order in the 
xy plane. Longitudinal order, present at A = 0, can be 
excluded in Fig. [T] for all U/t and for X/t > 0.002, so that 
the xyz AFMI phase is confined to a very small (possibly 
infinitesimal) interval starting at A = 0. The SM-TBI 
transition is evinced by the simultaneous opening of a 
single-particle and a spin gap, which as a function of A 
closely follow the U = results. The QSL-TBI transition 
for intermediate Hubbard U and small A turns out to be 
the most difficult and perhaps most interesting case, with 
the critical values extracted from a cusp (consistent with 
a closing) of the single-particle gap A sp and the spin gap 
A s . A more detailed discussion is given below. 

Our numerical results suggest that the TBI phase at 
finite U is adiabatically connected to the TBI state of the 



KM model (U = 0). Similarly, the QSL phase is stable 
over a finite range of A, in accordance with theoretical 
predictions.— Except for the smaller range of spin-orbit 
couplings compared to Ref. 0, which is chosen here to 
highlight the structure of the phase diagram around the 
QSL, we have obtained a number of additional points for 
the phase boundary of the QSL. The refined QSL phase 
boundary reveals a direct magnetic transition between 
the QSL and the AFMI phase at finite A. Our numerical 
data suggest the existence of a multicritical point where 
the QSL, TBI and AFMI phases meet. The estimated 
location of this point is (A c , U c ) rj (0.035i, 4.2i). 

Let us compare the phase diagram in Fig. [T] to other 
work. The magnetic phase boundary was calculated us- 
ing mean-field theory.— In that work, a transition from 
the TBI to an AFMI phase is observed, with the crit- 
ical U increasing with increasing A and comparable to 
the band width. However, the numerical values differ by 
up to a factor of two. The phase diagram from unbiased 
QMC simulations was presented by three of us»£ At that 
time, only one point on the QSL-TBI phase boundary 
was available, and the suggested dome-like structure of 
the QSL phase was based on the fact that the spin gap 
takes on its maximum around U/t = 4, in the middle of 
the A = QSL phase. Soon after this work, QMC results 
for the phase diagram were published by Zheng et al& 
Except for the absence of the QSL-TBI phase boundary, 
their phase diagram is compatible with previous^ and 
current results (Fig.Q}. The line X/t = 0.1 was studied by 
Yamaji and Imada using variational QMC simulations^ 
although with rather large quantitative differences con- 
cerning the location of the TBI- AFMI transition. The 
phase diagram has also been calculated using the varia- 
tional cluster approach^ and cluster dynamical mean- 
field theory»ii Apart from the fact that a true QSL phase 
is not accessible in any cluster calculation, the overall 
structure of the phase diagram in these works is con- 
sistent with Fig. [T] The quantitative phase boundaries 
seem to be slightly more accurate in the cluster dynam- 
ical mean-field case.— Both papers show a "QSL" -TBI 
phase boundary whose shape is in accordance with our 
refined phase diagram in Fig. [TJ 

Leei^ and Griset and Xui£ have recently made predic- 
tions about the nature of some of the phase transitions. 
In both works, the TBI-AFMI transition is argued to 
be in the 3D XY universality class, as already hinted at 
in Ref. 0. Griset and Xu further suggest that both the 
QSL- AFMI and the QSL-TBI transitions could be first 
order quantum phase transitions^ They also highlight 
the possibility of an additional, nematic order-disorder 
transition inside the AFMI at A = 0, instead of a pro- 
posed chiral AF order-disorder transition in the Hubbard 
model^ 9 . that should persist also at A > 0*^ The phase 
diagram of the KMH model has also been calculated us- 
ing analytical methods.—^ 

With the number of phases and their boundaries be- 
ing rather well established, the important open questions 
about the phase diagram concern the nature of the QSL 
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and AFMI phases, and of the various phase transitions. 
The structure of the phase diagram implies the existence 
of several distinct quantum phase transitions: SM-TBI, 
TBI-AFMI, QSL-AFMI, and QSL-TBI. The remaining 
SM-QSL transition only occurs at A = 0, and has been 
studied in detail before.— We discuss each of these tran- 
sitions below. 

The remainder of this section is organized as follows. 
We first consider the SM-TBI transition (at fixed U/t = 
2, see Fig. [TJ and the TBI-AFMI transition (at fixed 
X/t = 0.1), for which we can provide a fairly complete 
picture. From this we move on to the QSL-AFMI transi- 
tion (considering X/t = 0.0125), and finally the QSL-TBI 
transition (at U/t = 4). 



B. Semimetal to topological insulator transition 

We begin with the SM-TBI transition. To this end, 
we keep U/t = 2 fixed. In the absence of interactions, 
the spin-orbit term breaks the sublattice symmetry and 
generates a mass gap as well as a topological band struc- 
ture. Due to the underlying 17(1) spin symmetry, the 
band structure corresponds to two Haldane models with 
Chern numbers of opposite sign in the two spin sectors.— 

Figure [2ja) shows QMC results for the single-particle 
gap A sp and the spin gap A s as a function of A at U/t = 2. 
Starting in the SM phase at A = 0, where both gaps 
are zero, A sp and A s become nonzero for any finite A, 
and increase with increasing spin-orbit coupling. The 
QMC results at U/t — 2 closely follow the corresponding 
gaps for the noninteracting case U = 0,— A sp = 3-\/3A 
and A s = 2A sp . Interaction effects manifest themselves 
as a minor suppression of both gaps compared to their 
noninteracting values, especially at larger X/t, and by the 
spin gap falling below 2A sp . 

From these results, we draw the following conclusions. 
First, the SM phase of the Hubbard model is unstable at 
finite A, and hence only exists for A = 0, as indicated in 
Fig. [1] Second, the very small deviations in the depen- 
dence of the gaps on A compared to the noninteracting 
case suggest that the TBI phase at U > is essentially 
the same as at U = 0, provided U remains small enough 
to avoid the magnetic transition. This finding suggests 
that the two states are adiabatically connected. The mi- 
nor role of bulk interactions inside the TBI phase may 
be regarded as a consequence of the single-particle en- 
ergy gapfi and has been exploited to develop an effective 
model of the helical edges with a Hubbard U only at the 
edge sites of a ribboni^ 

The results for A sp and A s in Fig. HJa) are obtained 
from finite-size scaling, as shown for selected values of 
X/t in Figs.^Jb) and (c). Whereas A sp shows the famil- 
iar monotonic decrease with increasing system size, the 
spin gap reveals an unusual finite-size scaling behavior. 
Deep in the TBI phase [e.g., the top curve in Fig. [2Jc) 
corresponding to X/t = 0.03], A s systematically increases 
with increasing system size L. For small A (for example, 
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FIG. 2: (Color online) (a) Single-particle gap A sp , 2A sp and 
spin gap A s as a function of A at U/t = 2, across the SM- 
TBI transition. The QMC results are obtained from finite- 
size extrapolation using the fitting function Eq. (|21|) at A = 
and Eq. (|20p for A > 0. The numerical results shown as 
symbols agree well with the corresponding gaps of the KM 
model (U = 0)r A sp = 3^3 A and A s = 2A sp , revealing 
the adiabatic connection between the TBI at U = and at 
U > 0. (b) Finite-size scaling of the single-particle gap at 
selected values of X/t. (c) Finite-size scaling of the spin gap; 
for A > 0, we neglect the L = 3 results in the extrapolation. 



X/t = 0.0125) and small L, the scaling behavior is more 
complex, and only the system sizes beyond the crossover 
have been used in the extrapolation. The increase of A s 
with increasing L is a correlation effect; the spin gap is 
independent of L for [7 = 0. The observed increase of 
A s with system size can be reproduced using first-order 
perturbation theory in U. 

A possible physical explanation is inspired by the 
observation that for small L, A S (L) < 2A sp (L), see 
Figs. Hth) and (c). In contrast, the extrapolated val- 
ues almost match the relation for the noninteracting 



7 



case, A s w 2A sp , as shown in Fig. [2a). The strongly 
suppressed spin gap for small system sizes indicates 
pronounced particle-hole binding, driven by correlation- 
induced magnetic fluctuations. If the length scale of these 
fluctuations, which are a precursor of the magnetic tran- 
sition at U c , exceeds the system size, the spin gap is ex- 
pected to be suppressed similar to the magnetic phase 
where A s — > as L — > oo. Increasing L beyond the cor- 
relation length will restore the behavior expected for the 
weakly or noninteracting TBI phase. We will see below 
(Fig. [5]) that for larger U/t = 4, the spin gap is suppressed 
to values much below 2A sp even in the thermodynamic 
limit. Although a complete understanding of this effect is 
currently missing, we regard the unusual spin gap scaling 
as a signature of a correlated TBI. 

C. Topological insulator to antiferromagnet 
transition 

At large U/t, the TBI phase of the half-filled KMH 
model undergoes a transition to an AFMI, and TRI 
is spontaneously broken^ In the strong-coupling limit, 
U/t S> 1, the charge degrees of freedom are frozen and 
one can derive an effective spin model with antiferromag- 
netic nearest-neighbor Heiscnbcrg exchange J = 4t 2 /U 
that promotes isotropic magnetic order in the xy and z 
directions. The spin-orbit term of the KMH model re- 
duces the SU(2) spin symmetry of the Hubbard model 
to a U(l) symmetry corresponding to conservation of 
the total z-component of spin. Second-order perturba- 
tion theory gives an exchange interaction J' = 4X 2 /U 
between next-nearest neighbors. Importantly, the ex- 
change is antiferromagnctic in the longitudinal direc- 
tion, J'S?Sj, but ferromagnetic in the transverse direc- 
tion, —J'(SfSj + SVSj)& Combining all the exchange 
terms, magnetic order in the z direction becomes frus- 
trated, and the system favors an easy-plane Neel state. 
The so-called KM- Heiscnbcrg model was recently studied 
analytically^ 

From the above considerations, xy order is expected 
both for A = and A ^ 0. Hence, the phase boundary 
of the AFMI phase can be determined from the onset of 
transverse long-range magnetic order by monitoring the 
transverse structure factor 

a 

Kf1 q/3 = ^Z(-l) a (-l)^<*o|S+ S-,^ + S- a S+^ ) . 

rr' 

Here r,r' denote unit cells, a,j3£ {A, B} are sublatticc 
indices, (— 1) Q = 1 (—1) for a = A (B), and we have 
taken the trace of the corresponding 2x2 matrix of the 
structure factor. The quantity S^/N (with N = 2L 2 ) 
extrapolates to zero below U C (X), but takes on a finite 
value in the thermodynamic limit for U > U C (X), which 
is inside the AFMI phase of Fig. [T}£>£ It is also related to 
the transverse magnetization via m 2 — S^/N. 




1 /L 

FIG. 3: (Color online) Finite-size scaling of the rescaled mag- 
netic structure factor S^ F /N defined in Eq. (fTS")) at X/t = 0.1 
for different values of U/t, across the TBI-AFMI transition. 
The curves (representing polynomial fits) extrapolate to zero 
for U/t < 4.9, and to a finite value for U/t > 5.0, giving 
the critical value U c /t = 4.95(5). A more accurate estimate 
U c /t — 4.96(4) is obtained from Fig. [5ja). The inset shows 
the order parameter m xy as obtained from extrapolation to 
the thermodynamic limit. 



Numerical results for X/t = 0.1 are shown in Fig. [3] 
the transition is most obvious from the extrapolated or- 
der parameter shown in the inset. The extrapolation 
of Sj/^/N in system size gives a critical value U c /t = 
4.95(5). This value agrees with the slightly more accurate 
estimate U c /t — 4.96(4) which follows from the intersect 
of curves for different system sizes in Fig. [5ja). How- 
ever, this scaling analysis (see below for more details) 
relies on the knowledge of the universality class of the 
transition. By performing calculations at different X/t, 
we can determine the magnetic phase boundary, and we 
find good agreement with previous exact simulations at 
A = 0£ and A > OM. Our QMC results (not shown) fur- 
ther exclude the presence of longitudinal magnetic order 
along the entire TBI-AFMI phase boundary in Fig. Q] 
and up to U/t = 8. 

The TBI-AFMI transition is also reflected in the 
single-particle gap. Because the onset of long-range mag- 
netic order at the TBI-AFMI transition spontaneously 
breaks TRI, the transition from the TBI to the nonadia- 
batically connected AFMI can in principle occur without 
closing any excitation gaps. Instead, the transition man- 
ifests itself in A sp as a cusp at U c /t = 4.95(5), visible in 
Fig. 0Ja). The results are for the same value of X/t = 0.1 
considered in Fig. [3] A similar signature can be repro- 
duced already on the mean-field level, although with only 
a kink instead of a cusp at the critical point. Results for 
the gap and the mean-field order parameter are presented 
in the inset of Fig. HJa) . Figure Hta) also shows the clos- 
ing of the spin gap A s at U c ; the results were obtained 
from the finite-size scaling shown in the inset of Fig.[6][a). 

In Fig. 2][b), we present numerical data for the energy 
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FIG. 4: (Color online) (a) Single-particle gap A sp and spin 
gap A s as a function of U at X/t = 0.1. The values shown 
were obtained from an extrapolation to the thermodynamic 
limit. The dip in A sp and the closing of A s are consistent 
with Uc/t — 4.95(5). The inset shows the mean-field results 
for the single-particle gap and the magnetic order parameter, 
(b) Energy derivative with respect to U [Eq. ()16l) ] across the 
TBI-AFMI transition at X/t = 0.1 [U c /t = 4.95(5)]. 



derivative 



au [2 ^ 1 



(16) 



corresponding to the expectation value of the interaction 
term or, equivalently, the average double occupation, at 
X/t = 0.1. The continuous variation of this quantity 
across U c /t = 4.95(5) suggests a continuous transition. 

Having established the phase boundary of the magnetic 
transition at large U /t, we now consider the universality 
class. Given the remaining U(l) spin symmetry in the 
presence of spin-orbit coupling, the transition is expected 
to be in the 3D XY universality class. An intuitive pic- 
ture is based on local magnetic moments, which already 
exist in the magnetically disordered phase for U > 0, and 
order at U c . The onset of phase coherence at U — U c cor- 
responds to a 1/(1) symmetry breaking. This scenario is 
in accordance with the behavior of the spin gap A s in 
Fig. 2J The excitons are massive in the disordered phase 
(U < U c ), but condense in the ordered phase (U > U c ) 
where A s = 0. 

The conjectured 3D XY universality can be tested us- 




L U \U-U c )/U c 



FIG. 5: (Color online) Rescaled transverse magnetic struc- 
ture factor S^/N defined in Eq. (|15[) as a function of U at 
X/t = 0.1, for different lattice sizes L. Assuming the scal- 
ing form (T7J), (a) shows L 2f3/v S^/N . The intersection of 
curves for different system sizes yields U c /t = 4.96(4) for the 
critical point, (b) The scaling collapse obtained by plotting 
L W/u S ^i N ag a function of _ U c )/U c . The QMC 

data are fully consistent with the critical exponents 2=1, 
v = 0.6717(1) and /3 = 0.3486(1) of the 3D XY model.^ 



ing the zero-temperature, finite-size scaling forms 

SZ/N = L- 2 ?l v h[{U -UJL 1 '"] (17) 

and 



A s /t = L- z / 2 [(f/-f/c)i 1/iy ] 



(18) 



Here /i and /2 are dimensionless functions. The relevant 
critical exponents for the 3D XY model are z = 1, v = 
0.6717(1) and [3 = 0.3486(1)^ 

Using the same value X/t — 0.1 as before, we show in 
Fig. EJa) L 2 P/ U S X ^ F /N as a function of U for different 
system sizes L. If the scaling form Eq. (fTT)) with the 
critical exponents of the 3D XY model is correct, we 
expect to see an intersect of curves for different L at U = 
U c . As shown in Fig.[5{a), this prediction is indeed borne 
out by the QMC data, and we deduce U c /t = 4.96(4), 
in agreement with Fig. [3 Replotting L 2fi/v S^/N as 
a function of L X ^ V (\J — U c )/U c in Fig. EJb) produces a 
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v c ' c 

FIG. 6: (Color online) Spin gap A s as a function of U at 
X/t — 0.1, for different lattice sizes L. Given the scaling form 
Eq. (|18p . (a) shows L z A a . The intersection of curves for dif- 
ferent L gives U c /t = 4.96(4), consistent with Fig. [5ja). The 
inset shows the finite-size scaling of A s . (b) Scaling collapse 
obtained by plotting L z A s as a function of L 1// "((7 — U c )/U c . 
The QMC data are consistent with the 3D XY exponents 
« = l,i/ = 0.6717(1) and /3 = 0.3486(1)^ 
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FIG. 7: (Color online) (a) Finite-size scaling of the rescaled 
magnetic structure factor S^/N defined in Eq. (|15|) at 
X/t = 0.0125 for different values of U/t, across the QSL- 
AFMI transition. The data suggest a magnetic transition at 
Uc/t — 4.3(2). Lines correspond to polynomial fits, (b) Same 
as in (a) but showing the longitudinal structure factor S^f 
defined in Eq. (|19|) . In contrast to (a), there is no long-range 
order over the range of U/t values considered. 



D. Spin liquid to antiferromagnet transition 



clean scaling collapse onto a single curve. Figure [5] hence 
demonstrates that the assumption of 3D XY behavior is 
fully consistent with the QMC data. 

Figure IH1 shows a similar analysis for the spin gap A s , 
using the scaling form (|18|) . Although the statistical qual- 
ity of the data is not quite as good as for the structure 
factor, we again find satisfactory scaling (in particular, 
there is no noticeable drift of the intersect with increasing 
L) and the same U c using the 3D XY critical exponents. 

Based on the existence of a U(l) spin symmetry 
throughout the TBI phase, we expect the 3D XY behav- 
ior found at X/t = 0.1 to be generic for this transition, 
in agreement with previous predictions^ ' ' 15 The finite- 
size corrections to the 3D XY scaling behavior become 
more pronounced on approaching the possible multicriti- 
cal point, as verified explicitly for X/t = 0.05. According 
to Griset and Xu^ the observed 3D XY behavior at 
A > suggests the absence of a chiral AF order-disorder 
transition inside the AFMI phase even for A = 0^ 



The refined phase boundary of the QSL phase shown in 
Fig.Q]cstablishes the existence of a QSL-AFMI transition 
at finite A, in addition to the A = transition studied 
before^ Since the present work is concerned with the 
KMH model, we only consider finite values A > here. 

The simple picture of the magnetic transition as an or- 
dering transition of magnetic moments (or exciton con- 
densation) discussed in the context of the TBI-AFMI 
transition cannot straightforwardly be applied to the 
QSL-AFMI transition. For example, a Zi spin liquid ex- 
hibits charge fractionalization, and therefore has no well- 
defined magnetic modes. Fractionalization could lead to 
an unusually large anomalous dimension.— On the other 
hand, if the QSL phase was adiabatically connected to a 
simple band insulator (without charge fractionalization), 
the transition is again expected to be of the 3D XY type, 
similar to the TBI-AFMI transition. 

Due to the small size of the spin gap in the QSL phase, 
the extrapolation of the order parameter (ITS]) to the ther- 
modynamic limit is much more delicate than for the TBI- 
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AFMI transition. In particular, a scaling analysis along 
the lines of Figs. [5] and [S] is not conclusive with the cur- 
rently available system sizes. 

We first address the question of longitudinal magnetic 
order. The phase diagram presented by Yu et al.^ based 
on results from the variational cluster approach, shows 
an extended region inside the AFMI phase in which the 
authors claim that magnetic order exists both in the xy 
plane and in the z direction. For A = 0, this region is 
argued to extend all the way to the QSL-AFMI phase 
boundary, leading to a simultaneous onset of transverse 
and longitudinal order at U c . At A > 0, Yu et al*£ find 
a transition from the TBI to an xy ordered AFMI at U c , 
and an onset of z order at even larger values of U. Hence, 
for A > 0, there would be an additional crossover (no 
symmetry breaking) inside the AFMI phase. Whereas z 
order is known to exist in the Hubbard model (A = 0), 
this result is surprising in the light of the strong-coupling 
picture mentioned above, in which antifcrromagnctic cor- 
relations in the z direction are frustrated by the interplay 
of hopping t and spin-orbit coupling A. 

To clarify the situation, wc use unbiased QMC sim- 
ulations and calculate the transverse structure factor 
[Eq. (|15p] as well as the longitudinal structure factor 
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FIG. 8: (Color online) Energy derivative with respect to U 
[Eq. {TfiJ] across the QSL-AFMI transition at X/t = 0.0125 
[U c /t = 4.3(2)], and at X/t = 0.04 [U c /t = 4.3(2)], close to 
the multicritical point. There are no signs of a first-order 
transition. 



ra F r = 7iE(- i )"(- i ) /3 ^i^ 



at X/t = 0.0125. The results are shown in Fig.[7J The on- 
set of transverse magnetic order is visible from the finite- 
size extrapolation of S^/N depicted in Fig. [7|a), and 
the critical value U c /t = 4.3(2) is shown in the phase 
diagram in Fig. [TJ However, as revealed by Fig. Efb), 
there is no long-range order in the longitudinal direction 
even for large values of U/t = 8. We have carried out 
simulations down to X/t = 0.002, where longitudinal or- 
der would be most favorable, but found no z order for 
the U range shown in Fig. [TJ Hence, an extended re- 
gion of z order as suggested by Ref. [l0| does not exist, 
and the phase diagram is instead given by Fig. [TJ with a 
very narrow, possibly infinitesimal, region of coexisting 
longitudinal and transverse order near A = 0. The dis- 
crepancy between our exact numerical results and those 
of the variational cluster approach is most likely a conse- 
quence of the very small cluster sizes used for the latter. 
Although the strong-coupling picture with exchange con- 
stants J, J' is not justified for intermediate U, the frustra- 
tion in the z direction qualitatively explains the absence 
of longitudinal order found numerically. The purely in- 
plane magnetic order agrees with field-theory predictions 
for the KMH models 

Using field theory arguments, Griset and Xv*i£ sug- 
gested the possibility that the QSL-AFMI transition 
could be first order. To test this hypothesis, we show 
in Fig. [SJa) the energy derivative dF/dU [Eq. (|T5|) ]. Wc 
do not find any sign of discontinuous behavior near U c , 



which suggests that the transition is continuous. How- 
ever, we cannot exclude the possibility of a weakly first- 
order transition. We have also calculated dF/dU at 
X/t = 0.04 (close to the multicritical point) and found 
no signature of discontinuous behavior, see Fig. [8|b). 



E. Spin liquid to topological insulator transition 

A characteristic feature of both the QSL and the TBI 
phase is the absence of broken symmetries. Therefore, 
the QSL-TBI transition cannot be tracked by a local or- 
der parameter. In previous work)£ the critical point A c at 
U/t = 4 was determined from the behavior of the single- 
particle gap. Here we discuss the underlying procedure 
in detail, present new results with improved resolution of 
the critical point and based on larger system sizes up to 
L = 18, and refine the phase boundary by determining 
the critical point at two other values of U/t. Moreover, 
we show numerical results for the spin gap, and address 
the possibility of a first-order transition. 

FigurclHIa) shows the single-particle gap A sp as a func- 
tion of A at U/t = 4. The data points are obtained from 
extrapolation to the thermodynamic limit, as illustrated 
in the inset. Deep in the TBI phase (A > A c ), we use the 
fitting function (a = sp, s) 



A a (L)/t = a + e- L 't" {b/L + c/L 2 ) , 



(20) 



with a correlation length £ Q . On approaching A c from 
above, the correlation length set by the single-particle 
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FIG. 9: (Color online) (a) Single-particle gap A sp and (b) spin 
gap A s and inverse spin correlation length l/£' s as a function 
of A at U/t = 4, across the QSL-TBI transition. A sp is ob- 
tained from finite-size scaling using Eq. Q20p for X/t > 0.045 
and Eq. pi) for X/t < 0.045. The cusp defines the critical 
coupling X c /t — 0.030(1). The inset in (a) shows the scaling 
for selected values of X/t. A s is obtained from finite-size scal- 
ing using Eq. (|20|l . see (c). £ B is extracted from fits to the 
spin-spin correlation function defined in Eq. ()22|) at r — L/2, 
see inset in (b). 



gap, £ sp , increases and exceeds L/2 = 9 (L = 18 being 
our largest system size) for X/t w 0.04. For X/t smaller 
than 0.04, we use 



A a (L)/t = a + b/L + c/L 2 



(21) 



As a function of A, the extrapolated single-particle gap in 
Fig. inja) initially decreases when starting from the QSL 
at A = 0, reveals a cusp centered at X c /t = 0.030(1), and 
increases rather quickly with increasing A for A > A c . As 



in previous work^ we take the location of the cusp to 
define the critical point A c of the QSL-TBI transition. 
We will argue below that the data are consistent with a 
closing of the gap at A c , and that the cusp is a result of 
finite-size effects. For X/t > 0.045, the larger system sizes 
now available result in larger values of A sp compared to 
previous workil 

A similar analysis can be carried out for the spin gap 
A s using Eq. ([2U]) for A > A c . Similar to the SM-TBI 
transition discussed above, the spin gap shows an un- 
usual finite-size scaling inside the TBI phase, as shown 
in Fig.[SJc). The system sizes required to see saturation 
(i.e., the magnetic correlation lengths) are significantly 
larger at U/t = 4 than at U/t = 2, cf. Fig. H^c). The ex- 
trapolated values of A s for A > A c are shown in Fig.^b). 
Comparing A sp and A s [Figs. [9(a) and^b)], we see that 
in contrast to U/t = 2 [Fig. [2ja)] we have A s < 2A sp 
in the TBI phase at U/t = 4. The suppressed spin gap 
indicates substantial particle-hole binding. In the QSL 
phase, the small values of the spin gap make an accurate 
determination very challenging; A s is largest at A = 0, 
where it was previously determined as A s /t = 0.023(5) £ 

Figure H^b) reveals that the behavior of the spin gap 
for A > A c is very similar to that of A sp . As a consistency 
check, we also show the inverse spin correlation length £ 8 . 
Assuming the form 

S xx (r) = (S?SS) =e- r/ ^(a/r + b/r 2 ) (22) 

for the real-space transverse spin-spin correlation func- 
tion, and taking the largest available distance r = L/2 
for each system size [see inset of Fig. E^b)], the depen- 
dence of on A is in good agreement with A s and A sp . 
We only show the values ^ < L/2. 

For a noninteracting Z2 TBI, there is a simple rela- 
tion between the excitation gaps in the single-particle 
sector and, e.g., in the spin and particle- hole channels. 
In the presence of (strong) interactions, these relations 
may be modified, and we indeed find A s < 2A sp in the 
TBI phase above the QSL-TBI transition, as well as in 
the QSL phase at A = 0£ The argument that A sp has 
to close across a transition that involves a change of the 
topological indcxi holds only for the noninteracting case. 
In general, it is not clear which excitation gaps (one or 
more) close if the states on either side of the transition are 
not adiabatically connected. For example, in the inter- 
acting Haldane model^ there is an exact degeneracy of 
the three lowest states at the TBI to charge density wave 
transition. As a result, the first and second excitation 
gaps (E\ — Eq and E2 — Eq) close, but the single-particle 
gap A sp shows only a cusp at the critical point. 

As argued in previous work)! the results for A sp are 
consistent with a vanishing of the single-particle gap at 
A c . Furthermore, the results in Fig.[5]reveal that A sp and 
A s behave very similar on approaching A c from above, 
and we may therefore expect to see a simultaneous clos- 
ing of A sp and A s . Such a gap closing suggests differ- 
ent Chern numbers for the TBI and QSL phases. Ad- 
ditionally, the quick, almost linear opening of the gaps 
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for A > A c is reminiscent of Fig. [2ja) for the SM-TBI 
transition, suggesting that a non-TBI phase (the QSL) 
exists at small values of A, and that a transition to the 
TBI phase takes place at A c > 0. This picture confirms 
the expectation that the fully gapped QSL phase should 
be stable under a small perturbation in the form of the 
spin-orbit term. 

We attribute the small but nonzero values of the gaps 
at A c to finite-size effects. Although we used the same 
range of system sizes (up to L = 18) as for A = 0, 6 
the larger correlation lengths in the present case, espe- 
cially in the spin channel, make the analysis significantly 
harder. On approaching A c from above, the correlation 
lengths exceed the largest distance available on the clus- 
ters used. If we consider only the data points for which 
the correlation lengths fit on the largest system, i.e. the 
range X/t > 0.045 in Figs. [SJa) and^b), the functional 
form of A sp , A s and 1/^ strongly suggests a closing of 
the gaps very close to X c /t = 0.030(1). The fact that 
the finite-size scaled gaps in Figs. HJa) and [Hb) satu- 
rate at a finite value close to A c is therefore likely to be 
a result of insufficiently large system sizes. The latter 
require that we switch to the polynomial fitting func- 
tion ([21]) close to A c for A sp , and do not permit a reliable 
calculation of A s or ^ close to A c . The question if the 
gaps close or not cannot be answered using approximate 
cluster calculations ; 10 ! 11 because such methods are not 
capable of describing a true QSL phase. 

To determine the shape of the QSL phase boundary, 
we have calculated the single-particle gap for two other 
values of U/t. The critical values A c are again defined by 
the location of the cusp in A sp . We find X c /t = 0.025(2) 
for U/t = 3.8 and X c /t = 0.032(2) for U/t = 4.1. 

Recent theoretical work based on a 1/N expan- 
sion predicts the possibility of a first-order QSL-TBI 
transition^ To test this prediction, we show in Fig. [TU] 
the quantity 



dF 
dA 



(23) 



corresponding to the expectation value of the spin-orbit 
term in Eq. ([T]). For the range of system sizes, and on 
the very fine grid of A values, there is no sign of a dis- 
continuity. Again, we cannot rule out the possibility of a 
weakly first-order transition. 



V. CONCLUSIONS AND OUTLOOK 

Using exact quantum Monte Carlo simulations, we 
have obtained the phase diagram of the Kane-Mclc- 
Hubbard model (Fig. [IJ . For weak Hubbard interaction, 
the system is either a semimetal (SM) (at zero spin-orbit 
coupling, A = 0) or a topological band insulator (TBI). 
The latter is adiabatically connected to the noninteract- 
ing groundstatc of the Kane-Mele model, as evinced by 
the almost identical dependence of the single-particle and 
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FIG. 10: (Color online) Energy derivative with respect to A 
[Eq. ([23j)], across the QSL-TBI transition at U/t = 4. There 
is no sign of a discontinuity at \ c /t = 0.030(1). 



spin gaps on A. We have presented evidence for sub- 
stantial particle-hole binding in the TBI phase for small 
systems or large Hubbard interaction. For intermediate 
Hubbard U, the model supports a quantum spin liquid 
(QSL) phase at small A and a TBI phase at large A. At 
large U, long-range magnetic order breaks time reversal 
invariance, and the system becomes an antiferromagnetic 
Mott insulator (AFMI). In the presence of spin-orbit cou- 
pling, magnetic order is restricted to the xy plane. 

As previously suggested , 7 ' 14 ! 15 the magnetic TBI- 
AFMI transition can be understood as a condensation 
of magnetic excitons. A scaling analysis of the magneti- 
zation and the spin gap provides clear evidence for the 
3D XY nature of the transition. The onset of long-range 
order coincides with the closing of the spin gap, whereas 
the single-particle gap stays finite but shows a cusp at 
the critical point. In contrast to theoretical predictions, 
the corresponding transition between the QSL and the 
AFMI appears to be continuous. 

The QSL-TBI transition manifests itself as a cusp in 
the single-particle and spin gap. The numerical data are 
compatible with a complete closing of the gaps, but a def- 
inite conclusion is complicated by restrictions in lattice 
sizes. The independently deduced inverse spin correla- 
tion length is consistent with this picture, thereby sug- 
gesting that the QSL and TBI phases are not adiabati- 
cally connected. Finally, we find no sign of a predicted 
first-order transition. 

There remain a number of interesting open issues, in- 
cluding a characterization of the QSL phase, resolving 
the possible closing of the spin and single-particle gap 
across the QSL-TBI transition, and the universality class 
of the QSL-AFMI transitions both at A = and A > 0. 
Understanding the universality would provide important 
insight about the nature of the QSL phase, including the 
possible existence of fractionalization. All these ques- 
tions require significantly larger system sizes and hence 
massively parallel computers and will be addressed in fu- 
ture work. 
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